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Abstract 

The ionic atmospheres around nucleic acids play important roles in biological function. Large-scale explicit solvent sim¬ 
ulations coupled to experimental assays such as anomalous small-angle X-ray scattering (ASAXS) can provide important 
insights into the structure and energetics of such atmospheres but are time- and resource-intensive. In this paper, we use clas¬ 
sical density functional theory (cDFT) to explore the balance between ion-DNA, ion-water, and ion-ion interactions in ionic 
atmospheres of RbCl, SrCh, and CoHexCla (cobalt hexammine chloride) around a B-form DNA molecule. The accuracy 
of the cDFT calculations was assessed by comparison between simulated and experimental ASAXS curves, demonstrating 
that an accurate model should take into account ion-ion correlation and ion hydration forces, DNA topology, and the dis¬ 
crete distribution of charges on the DNA backbone . As expected, these calculations revealed significant differences between 
monovalent, divalent, and trivalent cation distributions around DNA. About half of the DNA-bound Rb+ ions penetrate into 
the minor groove of the DNA and half adsorb on the DNA backbone . The fraction of cations in the minor groove decreases 
for the larger Sr^+ ions and becomes zero for CoHex^+ ions, which all adsorb on the DNA backbone . The distribution 
of CoHex^+ ions is mainly determined by Coulomb and steric interactions, while ion-correlation forces play a central role 
in the monovalent Rb+ distribution and a combination of ion-correlation and hydration forces affect the Sr^+ distribution 
around DNA. This does not imply that correlations in CoHex solutions are weaker or stronger than for other ions. Steric 
inaccessibility of the grooves to large CoHex ions leads to their binding at the DNA surface. In this binding mode, first-order 
electrostatic interactions (Coulomb) dominate the overall binding energy as evidenced by low sensitivity of ionic distribution 
to the presence or absence of second-order electrostatic correlation interactions. 

Please address correspondence to Nathan Baker (nathan.baker@pnnl.gov). 


Introduction 

Interactions with ions stabilize nucleic acid secondary and tertiary structure, have a major impact on DNA packing in cells, 
and strongly influence protein and drug binding dD®. A fraction of counterions bind to specific sites on nucleic acids and can 
be detected in crystallographic structures (O, while other counterions form a dynamic ion atmosphere around DNA, diffus¬ 
ing along the molecule and exchanging with ions in bulk solution (HB- Mean-field approaches such as Manning counterion 
condensation (fTSl l and Poisson-Boltzmann (PB) (fTsHThl l theory have been used to obtain insight into ion distributions around 
biomolecules and ion-mediated interactions between macro-ions and have been compared with experimental data with some 
success (HMS. While successful in describing some properties of nucleic acids in electrolyte solutions (e.g., RNA pKa 
shifts (l2Tll . monovalent ion concentration linkages to ligand-DNA binding |9l |22l |2B . and low valency ion distributions 
around DNA), these mean-field methods often fail when the ion charge concentration increases. For example, PB models 
cannot capture the displacement of Na"*" by Mg^+ around DNA in mixed solutions (flTl l or ion-mediated DNA-DNA attractive 
interactions dSll. By imposing the constraint that a fraction of the counterions are bound (condensed) to polyelectrolyte and 
part form ionic atmosphere in the mean-field counterion condensation theory, it has been possible to reproduce attraction 
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Figure 1: Macroion models used in classical DFT simulations: (a) Model of a cylinder with uniform axial charge density; (b) 
discrete charge model. 


between like-charged polyelectrolytes in the presence of monovalent counterions in the intermediate range of separations 
(l24l428l l. Manning suggests that the origin of this effect lies in the increase in entropy due to the increase in the effective 
volume available for condensed counterions as two DNA molecules approach (l26l l. Such condensation implies penetration of 
ions through the DNA hydration layer and their partial desolvation to form direct bonds with DNA (fT8ll29llMb . Describing 
this process requires atomistic or coarse-grained representation of the macro-ion, which captures both the discreteness of 
charge distribution on the DNA backbone and DNA topology, as well as a model for ion desolvation. Such characteristics are 
not currently present in the PB equation or other popular models of biomolecular electrostatics. 

These failures suggest that, to reliably describe ion distribution around nucleic acids, the theoretical model must be rehned 
to include more detailed interactions and incorporate higher-order non-mean-held interactions such as Huctuations. Such 
extensions of PB approach have been developed for simple geometries (e.g., plates, rods, spheres, etc.) to include second-order 
terms representing the interactions between huctuations in ionic densities (l3TI4J7l) . These extended models and molecular 
simulations (l38ti46l l as well as experimental data (147145^ predict attraction between like-charged objects in the presence of 
multivalent electrolytes. 

In this study, we establish a minimal model based on classical density functional theory (cDFT) to systematically study 
the inhuence of the discrete DNA molecular charge representation, ion-ion correlations, and ion-solvent interactions on the 
distribution of monovalent and multivalent ions around highly charged macromolecules. We show that this model is able to 
accurately reproduce the results of anomalous small-angle X-ray scattering (ASAXS) experiments (119115914611) for B-DNA in 
RbCl, SrCl 2 and CoHexCla solutions. As expected, ion-ion correlations play a signihcant role in the accurate prediction of 
ASAXS curves. However, our results also demonstrate the importance of ion solvation in cation-DNA interactions and show 
that for doubly-charged cations these interactions can be as important as ion-ion correlations in modeling ion distributions 
around DNA. 


Methods 


DNA models 

We used two coarse-grained models for the DNA macro-ion in the cDFT simulations: an inhnitely long cylinder with a uni¬ 
form line charge density along its z-axis (charge distribution -1 e per 0.17 nm and the 2 nm cylinder diameter) and a model 
with a discrete charge distribution (Fig. [^. Na+ counterions present at 0.78 M concentration in all DNA calculations. The 
discrete charge distribution of the second model is described by three particle types: two helical arrays of charged spheres 
that represent the phosphate groups (charge -1 e, diameter 0.42 nm), two helical arrays of neutral spheres (diameter 0.42 nm) 
that represent the sugar/base groups, and an array of overlapping neutral spheres (diameter 0.78 nm) dehning the DNA axis 
([62). The positions of these spheres were chosen to mimic B-form DNA using a cylindrical coordinate system (r®, <pp^j) 
for DNA backbone s and base pair j. The phosphate spheres have coordinates r® = 0.89 nm, (j)j=(j)Q + 36j degrees, and 
Zj = Zq + 0.34j nm; the sugar/base spheres have coordinates r® = 0.59 nm, t/)® = -(- 36j degrees, and Zj = Zq + 0.34j 

nm; and the axis spheres have coordinates r = 0 nm, cj) = 0 degrees, and Zj = 0.5 4- 0.34j nm. There are 10 base pairs 
{j = 0,..., 9) per turn of B-DNA; the angular cylindrical coordinates for backbone start at = 0 and = 154 degrees, 
respectively. 
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Table 1: A summary of the different computational models used in this paper to assess the influence of different energetic con¬ 
tributions (ion-ion electrostatic correlations, ion-ion steric correlations, ion-solvent interactions, and water structural changes) 
on DNA-ion distributions and compare the resulting distribution functions with experimental ASAXS data. The rows provide 
model descriptions while the columns indicate which physical phenomena are included by the models. Additional quantitative 
dat a on model results are included in Table _ 


Model 

Ion-ion 

elec. 

correl. 

Ion-ion 

steric 

correl. 

Ion- 

solvent 

interact. 

Water 

struct. 

change 

Rb+ exp. 
agreement 

Sr^+ exp. 
agreement 

CoHex3+ 

exp. 

agreement 

NLPB 

no 

no 

no 

no 

no* 

no* 

yes 

cDFT, no correla- 

no 

yes 

no 

yes 

no* 

no* 

yes 

tion 

cDFT, no ion sol- 

yes 

yes 

no 

yes 

yes 

no* 

yes 

vation 

cDFT, full model 

yes 

yes 

yes 

yes 

yes 

yes 

yes 


*Agreement with experiment can be obtained by fitting ion radii. 


Computational models 

A variety of computational models were used with the DNA models described above to assess the influence of different ener¬ 
getic contributions on DNA-ion interactions. These models are summarized in Table[T]and described in detail in the following 
sections. 

Classical density functional theory (cDFT) 

Classical DFT (cDFT) has been previously used to determine the equilibrium distributions of multicomponent salt solutions 
surrounding DNA (l6Tl - l65l l. In our cDFT models, the aqueous salt solution was modeled as a dielectric medium with e = 78.5, 
charged spherical particles representing ions, and neutral spherical particles representing water molecules. The concentration 
of spherical “water molecules” was 55.5 M, chosen to model experimental water density. The solutions considered in this 
work were aqueous NaCl, RbCl, SrCl 2 and CoHexCla electrolytes in Na“*' buffer. These electrolytes were chosen based on the 
availability of the experimental data for these systems (I59ll60b . We used experimental crystalline ionic diameters for mobile 
ions: CTNa = 0.204 nm, crcoHex = 1.166 nm, ctsi- = 0.252 nm, cjRb = 0.322 nm, era = 0.362 nm, and CTwater = 0.275 nm 
(|66] |. The ion charges were gNa = +1, 9CoHex = +3, qsr = +2, (/Rb = +1, qci = -1, and gwatei- = 0. Parameterization of the 
cDFT model was performed against experimental data for the concentration dependence of mean activity coefficients in bulk 
electrolyte solutions (see Supporting Information). All calculations were performed at 298 K temperature. 

To determine the equilibrium water and ion distributions via cDFT, the total Helmholtz free energy functional is minimized 
with respect to the densities of all the species in the presence of rigid DNA models. For this optimization, it is convenient 
to partition the total free energy of the system into so-called ideal (J^“^) and excess components (J^“) (l63l l. The ideal free 
energy corresponds to the non-interacting system and is determined by the configurational entropy contributions from water 
and small ions, 

= {pfv)\ogpfv) - pfv))dv (1) 

I -In 

where k is Boltzmann’s constant, T is the temperature, pi -.Cl^ [0,1] is the density profile of ion species i, N is the number 
of ion species, r G fl is the ion coordinate, and H G is the calculation domain. The excess free energy is generally not 
known exactly but can be approximated by 


j-ex 


TO I -pi 

hs ' f 


TT^x 

corr '' C 


+ P. 


solv 


( 2 ) 


where is the hard-sphere repulsion term, is the ion-ion electrostatic correlation term, is the direct Coulomb term, 
and is the ion-solvent interaction term included in some cDFT calculations (as described below). 

The ion-ion steric correlation term excess free energy describes ion and water many-body interactions in condensed phase 
due to density fluctuations and can be approximated by Fundamental Measure Theory (l67l) as 




hs 


[rin, (r)] dv 


(3) 
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where the functional has the form (l68T l 


^’hs(r') = -no ln(l - 713) + + 

1-713 


8677773 


ln(l - 773) + 


ni • 712 


8677773(1 - n3)2_ 


772 


1-773 


1277 Tl^ 


ln(l - 773) + 


1277773(1 - 773)2 


772(712-712). ( 4 ) 


where Ua and tip are the scalar and vector weighted averages of the density distribution functions pi{r) and are defined by: 

'na{r)='^[ pi{r')uj\°'\r'- r)Ar', fora = 0,l,2,8 

i 

Pi{r')i^f\r'-r)dr', for/3 = l,2. 
r2 


In the limit of a bulk hard-sphere fluid in the absence of external fields, vector densities 711 and 712 vanish. In the same limit, 
the four scalar weighted densities reduce to the sum of bulk densities for all species ( tiq ) and the ID (777), 2D (772), and 3D (773) 
packing fractions. The “weight functions” and characterizing the geometry of particles (ion-ion steric correlations 
with radius Ri for ion species i), are given by (l6^ 


"w 

II 

Qi 

1 


(5) 

wf)(r) = |V6l(|r| - i?,)| = (5(|r| 

-R^) 

(6) 

= V6»(|r| - R^) = ^(5(|r 

-R^) 

(7) 

ujf\r) = a;fl(r)/(47ri?2) 


(8) 

w(i)(r) = a;fl(r)/(47ri?,) 


(9) 

a;fl(r) = wf l(r)/(47rii,). 


(10) 


In the preceding formula, 9 is the Heaviside step function, with 9{x) = 0 for a; > 0 and 9{x) = 1 for a; < 0, and S denotes 
the Dirac delta function. 

The ion-ion electrostatic interaction term can be derived using the Mean Spherical Approximation (l^lTSl l 


r ^ 
o 


kT 


r r ^ 

J Q. J Q. ^ j 


where is the bulk concentration of ion species i and the first term describes ion correlation free energy in bulk electrolyte 
solution in the absence of DNA. The first-order direct correlation functions are defined as 


(1) ^ k-i 
kT' 


( 12 ) 


where pi is the chemical potential of ion species i. The second-order direct correlation functions are defined as 



where qi is the charge of ion species i, e is the dielectric constant of the solvent, try = (<77 + cry) /2 is the hard-sphere contact 
distance between ions of diameters ai and cry, B is given by 
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^ = naij, K is the inverse Debye length = Ib = -gjq is the Bjerrum length, and e is the unit charge. The 

direct Coulomb free energy term can be calculated exactly 


-^c 



|r-r'| 


Pi{r)pj{v')dTdY'. 


Finally, the ion-solvent interaction term models ion-water interactions with a square well potential 


(15) 


V(r — r') 


{ CO |r-r'|<(Ty 

—£ aij < |r — r'l < -I- h 

0 CTy -I- /i < |r - r'l, 


(16) 


where e is the well depth, and h is the well width. For the current study, h = 0.2<; is the well width for interactions between 
ions and water and is the sum of radii of interacting particles dH. The following well depths were calculated using SPC/E 
water using the parameters from Horinek et al: £s,- = 0.01038 eV, sq = 0.0053894 eV, £Rb = ccoHex = 0.0021 eV (l69]l. Sim¬ 
ulations of concentration dependence of ion activity coefficients in RbCl and CoHexCla solutions demonstrated that adding 
attractive ion-water interactions does not affect the ion chemical potential (see Supporting Information). 

Minimization of the excess free energy functional with respect to the water and ion densities gives 


Pi(r) = exp 



1 \ 
kT 5 p^{y)) ■ 


We solve Poisson’s equation 

- V • e(r)V(^(r) = ^ q^Pi{v) 

i 


(17) 


(18) 


for the electrostatic potential (y3(r)) where e(r) is the dielectric coefficient. For an infinitely long uniformly charged cylinder 
in electroneutral conditions, the potential 


AtT /T\ ^ 

^(r) = ^J t\og(J^jY^q,p^{t)dt. 


(19) 


Using this potential for the cylinder model and a numerical solution to Poisson’s equation (Eq. 181 for the 3D DNA model, 
the expression for the densities is 


p,(r) = exp ( g; - 


qMr) 1 + + 


kT 


kT 




( 20 ) 


The resulting system of Eqs. 18 and 20 was solved iteratively to self-consistency using the numerical procedure described in 
detail by Meng ( |68] |. In particular, equilibrium ion density distributions were obtained using a relaxed Gummel iterative pro¬ 
cedure for 3D systems and Picard iterations in ID. Convergence was considered to be achieved when the maximum difference 
between the input and the output density profiles between iterations was smaller than 10“®. The solution of Eqs. 18 


and 


20 


encompasses the equilibrium distribution of the densities of all ion species, corresponding to the minimum of the total free 
energy; the corresponding free energies for each contribution; and the chemical potentials. “Panoramic” density distributions 
representing angular distributions of ions on DNA backbone and in minor grooves were calculated along the cotTesponding 
helical shells. Eor each, angle the ion densities were averaged within the shells over r and z. Eor ions on the DNA backbone 
and in the minor grooves, the radial positions of the shells was defined as 1 < r < (1 -f cr) nm and 0.5 < r < 1 nm, 
respectively. 

Three main features distinguish our approach from previous cDET models (l70l472l l. Eirst, our model includes a full rep¬ 
resentation of the coarse-grained DNA topology and a discrete distribution of charges. Second, we use Pauling diameters for 
ions and van der Waals diameters for water molecules as opposed to previous restricted models where all species have the 
same diameter. Einally, our model includes water-ion attractive interactions. 
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Table 2: Numbers of condensed ions for the models described in Table [T] and experimental results (l59l l (where available). 


Model 

Rb+ 

Sr2+ 

CoHex3+ 

NLPB 

23.0 

13.6 

5.6 

cDFT, no correlation 

23.4 

13.7 

5.6 

cDFT, no ion solvation 

34.9 

25.1 

5.7 

cDFT, full model 

34.6 

16.6 

5.7 

Experiment 

34 ±3 

19 ±2 

No data 


Anomalous small-angle X-ray scattering curve calculations 

ASAXS profiles were calculated using the ion density distributions pi{r) around DNA. In the 3D model, ion densities were 
averaged in cylindrical coordinates over the cylinder azimuthal angle (f) and length z for each radial distance r from the DNA 
axis. The excess form factor for ion species a was calculated as 

Fion,aiQ) = J Pair)e~'^^^dr, (21) 

where Oq, is a constant related to the average electron density of ion species a and Q is the scattering vector. In the current 
study, we only consider the excess form factor due to cation species; the chloride anion has no ASAXS response. Furthermore, 
we only consider a single cation species at a time so that Flon(Q) = Fion,aiQ)- 

The excess form factor of DNA (FbNA(Q)) was calculated using AquaSAXS (l73l l. from the form factor of DNA in vacuo 
(^dna(Q))’ '^he form factor of the volume of water excluded by DNA (F’ona(Q))’ *^he form factor of hydration shell of 
the DNA (Fh,h(Q)): 

FunaIQ) = i^DNA(Q) - PwF^^iiQ) + p„Fhsh(Q), (22) 

where p^, is the bulk density of water. The form factor of the hydration shell is calculated using water density maps, 
Phsh(r), obtained via AquaSol (l74l) . which employs the Poisson-Boltzmann formalism with water treated as an assembly 
of self-oriented dipoles: 

Fhsh(Q) = b J e-^Q^dr (23) 

where 6 is a scale factor to adjust the hydration shell contribution (usually b = 1.0) and integration is performed over the 
region where solvent density deviates from the bulk by a factor larger (in magnitude) than ±10“'*. 

The ASAXS intensity is then calculated from these quantities as 

I{Q) = 2 ifUEl) - fLiE2)) {fDNAN,onFj,MQ)EUQ) + fionoNLFUQf) + (fZiE,) - fZiE^)) N^onELiQ) 

(24) 

where f[o„{Ei) is the energy-dependent real part of ion anomalous scattering factor,i?i is the energy far from the X-ray 
absorption edge of the ion, E 2 is the energy near the edge where ion scattering is suppressed by absorption, /jono is the energy 
independent solvent-corrected scattering factor, /dna is the effective number of electrons from DNA and is the number 
of excess ions around DNA (l59l l (see Supporting Information for more details). Since experimental data are available in arbi¬ 
trary units, theoretical intensities were uniformly scaled with a common scaling factor, chosen to match the experimental and 
calculated intensities, obtained using 3D cDFT-full model, at low Q. 


Results 

Tables [T] and [^provide qualitative and quantitative results on the performance of the models. These results are described in 
greater detail below. 


Comparison between DNA Model systems 

The uniformly charged cylinder model (Fig. [T] left) represents a one-dimensional case for which ionic distribution is only a 
function of the radial distance from the cylinder axis. This ID model produces monotonically decreasing with the distance 
from the cylinder surface density distributions of monovalent and multivalent counterions (Fig. |^. Competitive cation con¬ 
densation in mixed 5 mM CoHexCla - 1 - 20 mM NaCl solutions results in preferential CoHex^+ condensation on the cylinder 
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Figure 2: Ionic distributions around a uniformly charged cylinder. The distance x/d is the distance (x) from the cylinder 
surface scaled by the ion radius (d). (a) Solutions of 100 mM NaCl, 100 mM RbCl, and 10 mM SrCl 2 in a 1 mM NaCl buffer; 
concentration prohles are shown for Na+ (solid line), Rb+ (dashed line), and Sr^+ (dot-dashed line), (b) Solutions of 5 mM 
CoHexCla in a 20 mM NaCl buffer; concentration prohles are shown for Na+ (solid line), CoHex^+ (dashed line), and Cl“ 
(dot-dashed line). 


surface: sodium ions are not found in the immediate vicinity of the DNA (Fig. |^). This competition is in qualitative agree¬ 
ment with experimental observations of a negligible effect of Na+ on CoHex^+ binding when NaCl concentration is below 
40mM(|25]l. 

For monovalent ions, the ID cDFT calculations predict 91.5% and 77.5% DNA charge neutralization by Na+ and Rb+, 
respectively. Such differences in monovalent cation condensation on DNA were not observed experimentally (l7^ . demon¬ 
strating a fundamental dehciency of a uniformly charged cylinder model for simulating ionic atmosphere around DNA. For 
divalent ions, the ID cDFT calculations predict charge inversion at the DNA surface in SrCl 2 solution. Note that charge 
inversion in the presence of multivalent salts has also been observed in cDFT and MC simulations for a cylinder DNA model 
Go] ED- However, we do not see this effect in our more detailed 3D DNA geometry simulations (see below). Finally, for 
trivalent ions, 90% DNA charge neutralization is found within 5 CoHex radii from the cylinder surface or within the region 
where excess CoHex concentration is present (Fig.|^). 

We also performed 3D cDFT calculations of the same electrolyte solutions surrounding the helical discrete charge model 
(Fig. [fright). Fig. shows cDFT results for the monovalent ion Rb+. As shown in the panoramic density prohles, cDFT 
predicts a two-peak radial density distribution of Rb+ (Fig. [^: hrst peak at around 0.6 nm is due to cation penetration into 
DNA minor grooves and the second peak at 1.2 nm to Rb+ condensation on DNA backbones. These data are in good quanti¬ 
tative agreement with molecular dynamics results obtained using TIP3P water model (fTSl l. Fig. |^shows cDFT results for the 
divalent ion Sr^+. In the case of Sr^'*', the effect of ion-solvent interactions can be clearly seen in the density distribution of 
Sr^+ with respect to the DNA axis (Fig.|^. While both cDFT models - with and without ion-solvent interactions - produce 
two-peak Sr^+ density distributions at the same positions with respect to the DNA axis, the density distributions are qualita¬ 
tively different. In particular, the model without ion-solvent interactions predicts much higher Sr^+ concentration in the DNA 
grooves than on the backbone , while the trend is reversed in the model with ion-solvent interactions. The 3D cDFT results 
for trivalent CoHexCla solutions are shown in Fig.|^ and are very similar to those obtained from the ID cDFT model. 

Finally, we used the results of our cDFT and NLPB calculations to determine ASAXS prohles as described in the Meth¬ 
ods section. The results of these calculations for Rb’*’ and Sr^+, together with experimental data, are shown in Fig.|^ Similar 
results for CoHexCls are shown in Fig. 
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Figure 3: Cation distributions around DNA calculated using (a) cDFT with the full model, (b) cDFT with no ion-correlation 
interactions (cDFT-nc, dotted line), and (b,c) NLPB. In panels (a-c), panoramic Rb+ density distributions are shown on the 
DNA backbone (black lines) and in the minor groove (red lines) as defined in the manuscript text. The inset in panel (b) shows 
a zoom-in into a low-density region. Panoramic views of cation distributions around DNA in 100 mM NaCl are shown in blue 
in panels (a) and (c) for comparison. The radial Rb+ density distribution calculated from the full cDFT model is shown in 
panel (d). 


Discussion 

Comparison to Manning condensation 

The cDFT calculations of ionic distributions for the uniformly charged cylinder model (Fig. reproduce the Manning con¬ 
densation limits ([nil with approximately IM concentrations of singly-charged cations at the cylinder surface. A complete 3D 
cDFT model also reproduces the Manning condensation limit for monovalent cations: the concentration of condensed Rb+ 
on the backbone and in minor groove is about 1 M (Fig. [^. The 1 M limit for monovalent ion concentrations at DNA sur¬ 
face corresponds to 76% compensation of native B-form DNA charge by condensed counterions. Manning’s theory predicts 
that the concentration of condensed counterions is independent of bulk salt concentration in the range of 0.0001-0.1 M and 
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Figure 4: Panoramic density distributions of Sr^+ ions on (a) the DNA backbone and in (b) the DNA minor grooves obtained 
10 mM SrCl 2 . Backbone and minor groove definitions are provided in the text. Results are shown for the full cDFT model 
(solid red line), cDFT with no corrlations (cDFT-nc, dotted line), and the NLPB model (solid black line). Panel (c) shows 
radial densities of Sr^+ ions around a DNA molecule calculated using the full cDFT model (solid line) and the cDFT model 
without ion-water attractive interactions (dotted line). 


increases slightly for higher ionic strengths of monovalent electrolyte solution (80% charge compensation for 0.5 M and 83% 
for 1 M solutions). Additionally, the 3D cDFT model predicts that the multivalent ions form much denser layers at the DNA 
surface than the monovalent cations (Figures and|^, consistent with Manning theory. The good correlation between 
our non-mean-held cDFT model (with full ion-ion correlations included) and the mean-held Manning theory is somewhat 
surprising, particularly given the signihcant differences observed in the total condensed ion densities between cDFT and the 
mean-held NLPB approaches. However, Manning theory indirectly accounts for interactions beyond hrst-order electrostatics 
through partitioning the total ion density into condensed ions and the surrounding ionic atmosphere. This accounts for the 
success of Manning theory in predicting the condensed concentrations of 1:1 electrolyte counterions on DNA as observed in 
experiments (fT^ and recent MD simulations (fTSl) . 
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Figure 5; (a) Panoramic density distributions of CoHex^+ ions on DNA backbone obtained using the full cDFT (solid red 
line) model, the cDFT model with no correlations (cDFT-nc, dotted line), and the NLPB model (solid black line) for 0.5 
mM CoHexCla. The NLPB and cDFT-nc curves have been shifted by 0.5 mM for clarity, (b) Radial density distributions of 
CoHex^+ ions around the DNA molecule. Note: as seen in panel (b), there is zero CoHex density^"*" in the minor groove, so 
the corresponding panoramic density is not shown. 


Ion interaction with DNA grooves 

By dehnition, the cylinder model does not allow ion penetration inside DNA and therefore yields well-known monotonically 
decreasing counterion distributions shown in Fig. Thus, the model is not adequate for describing the interaction between 
DNA and small weakly solvated Na+ and Rb+ ions, which are known to penetrate into the minor grooves of DNA (fTSliTTl l. 
However, our more detailed helical charge model allows ion penetration. Simulations of RbCl solutions using this model 
showed that about half of the condensed Rb+ ions are bound to the minor groove of the DNA molecule (Fig. [^. The dis¬ 
tributions of cation densities on the DNA backbone and in the minor groove are highly structured: they exhibit a periodicity 
correlated with the periodic spacing of phosphate groups on the DNA backbone . In contrast, cation distributions in the major 
groove are mostly featureless (see Fig. S1 in Supporting Information), in agreement with previous simulations and experi¬ 
mental data (l30l l. Penetration of some cations into the grooves lowers the effective charge density on the DNA, limiting cation 
condensation on the backbone . 

Increasing cation valency correlates with a stronger preference of cation binding to phosphate groups on the DNA back¬ 
bone (Figures 1^1^ and[^. A similar preference for CoHex^+ binding to phosphates of B-DNA was also observed in MD 
simulations (TtM ^ and is determined by the strong electrostatic attraction of the trivalent cations to phosphate groups, CoHex^"''- 
CoHex^+ repulsion, and steric inaccessibility of B-DNA minor groove to the large CoHex^+ ions. As shown in Figuresj^and 
1^ both Sr^+ and CoHex^+ ions preferentially bind to every fourth phosphate on the strand. Further away from the DNA 
axis, the Sr^+ density variations along the angular cylindrical coordinate have the same period as the period of the angular 
phosphate distribution (Fig.|^. The period of the density variations for CoHex^+ (Fig.|^ is two times larger than for Sr^+; 
i.e., some Sr^+ions can penetrate into the minor groove, while CoHex^+ ions bind exclusively to phosphate groups on the 
backbone . 


Influence of correlation on ion distributions 

To investigate the influence of ion correlation forces on the distribution of ions around DNA, we used a cDFT model with¬ 
out ion-correlation interactions (cDFT-nc) as well as the nonlinear Poisson-Boltzmann (NLPB) model, which also lacks 
correlation (Figures and [^. Both models without correlations yield qualitatively different ion distributions than the 
3D cDFT calculations which include correlations. In the presence of correlations, sterically allowed ions accumulate in the 
minor groove; in the absence of correlations, ions accumulate near phosphate groups on the exterior of the DNA strand. The 
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Figure 6: Simulated and experimental ASAXS profiles for 25 bp DNA in (a) 100 mM RbCl and (b) 10 mM SrCl 2 solutions. 
Experimental data (I59ll60l) are shown as black dots. This figure shows simulation results using the full 3D cDFT model (blue 
lines), the 3D cDFT model without ion-water interactions (red lines), the full ID cDFT model (blue dots), and the NFPB 
model (green lines and dots). 



Figure 7: Simulated and experimental ASAXS profiles for 25 bp DNA in 0.5 mM CoHexCla DNA solutions. Experimental 
data (unpublished) are shown with a thick black line and those from Andresen et al dbOl l with a blue line. ID cDFT results are 
shown with a thin black line and ID NFPB results with a red line. The 3D cDFT and NFPB data coincide and are shown with 
a green line. 


largest qualitative difference between NFPB and cDFT ion distributions was observed for the Rb+ density distribution. In 
NFPB, Rb+ ions decorate the phosphate groups, driven by Coulombic interactions; the panoramic distribution of Rb+ ions 
condensed on the backbone in NFPB model has a larger peak at 45° and a smaller one at 135° (Figures and|^. In 
contrast, ion-ion correlations reduce the effective electrostatic repulsion between cations promoting their penetration into the 
grooves.Due to stronger Coulomb interactions between multiply-charged cations the effect of correlations is weaker for Sr^+ 
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and CoHex^+ resulting in the decrease in the fraction of counterions in the grooves with ion radius and charge (Fig.|^and Fig. 
I^. For CoHex^+ the concentration of counterions in the grooves becomes insignihcant.As a result, NLPB and cDFT predict 
qualitatively similar panoramic density distributions on DNA the DNA backbone for Sr^+ and CoHex^+ (Figures]^ and |^. 

As illustrated in the hgures and Tables [T] and the models without correlations (NLPB and cDFT-nc) are very simi¬ 
lar to each other, indicating the major influence of correlation on even low charge-density (monovalent) ion behavior. This 
result contrasts the conclusion that correlations are insignihcant in monovalent electrolytes from early theories of ion cor¬ 
relations ®. However, these theories considered electrolytes at uniformly charged surfaces ignoring the influence of the 
discreteness of charge distribution on fluctuations in ionic atmosphere. Not surprisingly, these models do not capture the exper¬ 
imentally observed attraction between like-charged polyelectrolytes in low concentration monovalent electrolytes l54b . 
Recent molecular dynamic simulations also point to the importance of non-mean-fleld interactions between biomolecules and 
monovalent electrolytes manifested in a more structured ionic atmosphere than that predicted by NLPB (fTSllTTb . The small 
difference between the cDFT-nc and NLPB models (in the height of the double peak around 45°) for Rb+ is due to the solvent 
excluded-volume effects included in the cDFT-nc model and absent from NLPB theory (Fig.|^. 

Ion correlations also influence ion-specific details in density distributions for counterions of the same valency. 3D cDFT 
results demonstrate that smaller Na'*' ions tend to accumulate on the DNA backbone and minor groove while Rb+ ions are 
more evenly distributed along the DNA helix (see Fig.[^). In the cDFT-nc and NLPB models with no correlation, the differ¬ 
ences between Na+ and Rb+ distributions are significantly smaller (see Fig. [^), suggesting that ion-correlation interactions 
are responsible for this effect. 


Comparison with ASAXS experiments 

ASAXS profiles calculated using the 3D cDFT model show very good agreement with experimental data (l60l l for RbCl solu¬ 
tions (Fig.[^. The shapes of the scattering curves are very similar in the ID cDFT and 3D cDFT models, with similar average 
numbers of condensed counterions; the ID and 3D cDFT calculations predict 34.9 and 34.6 condensed Rb+ ions, respectively. 
Both predictions are within error of the experimental measurement of 34 ± 3 ions (EIll. However, the distribution of Rb+ ions 
is different in these models: all condensed cations decorate the cylinder surface (by definition) in the ID model while half of 
the condensed cations are on the DNA backbone and the other half are in minor grooves in the 3D models. In contrast, the 3D 
NLPB model shows a significant deviation of the simulated scattering curve from the experimental one (Fig.|^. As discussed 
in the previous section, penetration of some cations into DNA grooves reduces the negative electrostatic potential acting on 
cations in solution. In the absence of any interactions beyond Coulomb forces, this penetration leads to lower concentrations 
of cations on the DNA surface and lower concentrations of condensed counterions. Previous NLPB simulations demonstrated 
that adjusting the ionic radius of Rb+ to its hydrated radius and prohibiting ion penetration into the DNA hydration shell can 
lead to closer agreement between 3D NLPB results and experiment (l76l l - but at the price of an incorrect ion distribution 
around DNA. 

Comparison of calculated and experimental SrCl 2 data highlight the importance of ion-solvent interactions on ion dis¬ 
tributions around DNA (Fig. |^. Sr^+ ions have a significant hydration energy; approximately 3 times higher than that of 
monovalent alkali metal ions. These strong cation-water interactions lower the entropy of water molecules around cations, 
but introduce a higher enthalpy cost for partial ion desolvation (ISOl l. Simulations with the solvent approximated as dielectric 
continuum do not account for such desolvation, limiting interactions in the system to first- and second-order electrostatic inter¬ 
actions: Coulomb and ion-correlation forces. To understand the importance of these desolvation contributions, we used two 
variants of the cDFT model; one with attractive cation-water interactions and another without. As shown in Fig. ASAXS 
curves calculated using cDFT without desolvation contributions via ion-solvent interactions deviate significantly from the 
experimental data and the ASAXS curves calculated using the complete cDFT model. On the other hand, the experimental 
ASAXS curves agree with those calculated from the complete cDFT model. The importance of desolvation is also empha¬ 
sized by the fact that inclusion of such interactions is essential for reproducing the chemical potentials of divalent cations 
but is not required for weakly hydrated alkali metal ions or CoHex^+ (see Supporting Information). In summary, ion-solvent 
interactions are important for accurately modeling ion-DNA interactions: desolvation reduces the excess chemical potential 
of cations and anions, lowering the effective concentration of electrolyte and weakening ion-ligand interactions. Surprisingly, 
an NLPB model that includes neither ion-solvent interactions nor ion-ion correlations reproduces the experimental ASAXS 
curves for Sr^+ (Fig. [^, although some differences are obvious in the more detailed radial density functions (Fig. I^. This 
agreement is serendipitous and is due to cancellation of errors from the lack of ion-correlation, that favors ion accumulation 
in the grooves, and ion solvation, that limits ion concentration in the grooves. 

Unexpectedly, the trivalent CoHexCla solution is the simplest ion to model around DNA; CoHex^+ can be reliably 
described by first-order electrostatics (i.e., direct Coulomb interactions). CoHex^+ ions decorate DNA backbone phosphates 
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and do not penetrate inside B-DNA grooves (Fig.|^. All models explored in this paper show reasonable agreement between 
the calculated ASAXS profiles and the experimental data (Fig.|^, and is consistent with all-atom MD simulations in explicit 
solvent (fTST l. For these triply-charged ions, the good agreement between cDFT and NLPB is caused by the dominance of 
first-order electrostatics in ion-DNA interactions over higher-order ion-ion correlations. Moreover, the large CoHex^+ ionic 
diameter creates a steric barrier for ion penetration inside the grooves, rendering the ID cylinder models adequate for cal¬ 
culating the average number of condensed CoHex^+ ions. Finally, because the diameter of CoHex^+ is large, the field at its 
surface is comparable to Na+ and the effects of ion-solvent interactions are lower than for the smaller divalent Sr^+ ions. It 
follows from the current study that the models required to describe CoHex^+ around a single B-DNA backbone are relatively 
simple. However, we expect that ion correlation forces will dominate DNA-DNA interactions between multiple backbones 
due to the higher local phosphate charge density. 

Conclusions 

We have studied the details of ionic atmospheres around DNA molecule for 1:1, 2:1 and 3:1 electrolytes using a combination 
of cDFT and NLPB methods. Our calculations demonstrated that ion-ion correlation interactions induce counterion penetra¬ 
tion into the DNA grooves, unless sterically prohibited by large ion radii. In particular, ion binding in the grooves - compared 
to binding on the backbone - has a profound effect on ion-induced nucleic acid condensation as demonstrated in our previous 
work (|7^ . Ion-solvent interactions have an opposite effect: when the enthalpy cost of desolvation is high (e.g., for Sr^+ ions), 
ion-water interactions limit ion penetration into the DNA grooves. Partial compensation of these two opposing effects explains 
the success of NLPB in reproducing the average number of condensed cations and the shape of the ASAXS curves of the ion¬ 
counting experiments. In contrast, cDFT model without ion-desolvation interactions was found to systematically overestimate 
ion concentration in DNA grooves. Through the comparison of several cDFT models and experimental data, we demonstrated 
that a minimum model to describe ion-polyelectrolyte interactions should include long-range correlations arising from density 
and charge density fluctuations in electrolyte solution as well as short-range ion (de)solvation forces. The latter interactions 
are often ignored in reduced models of electrolyte solutions limiting their applicability to the classes of weakly solvated ions. 
Ion hydration forces are particularly pronounced in solutions of multiply-charged ions and give significant contribution to 
ion activity and, therefore, to ion-polyelectrolyte interactions. Our results highlight important aspects of the properties of 
electrolyte solutions influencing ionic atmosphere around biomolecules that may significantly impact DNA condensation and 
biomolecules-ligand interactions. One caveat of the current work is its neglect of DNA sequence-specific effects which have 
been shown to influence ion binding in some cases dsBiHi. The goal of our initial research was to understand the general 
characteristics of DNA-ion interactions that drive the behavior of different ionic species around DNA. In the future, we plan 
to extend the DNA model to include such sequence-dependent structural variations. 
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3 Cation distributions around DNA calculated using (a) cDFT with the full model, (b) cDFT with no ion- 

correlation interactions (cDFT-nc, dotted line), and (b,c) NLPB. In panels (a-c), panoramic Rb'*' density 
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in the manuscript text. The inset in panel (b) shows a zoom-in into a low-density region. Panoramic views of 
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4 Panoramic density distributions of Sr^+ ions on (a) the DNA backbone and in (b) the DNA minor grooves 
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5 (a) Panoramic density distributions of CoHex^+ ions on DNA backbone obtained using the full cDFT (solid 
red line) model, the cDFT model with no correlations (cDFT-nc, dotted line), and the NLPB model (solid 
black line) for 0.5 mM CoHexCla. The NLPB and cDFT-nc curves have been shifted by 0.5 mM for clarity. 

(b) Radial density distributions of CoHex^+ ions around the DNA molecule. Note: as seen in panel (b), there 

is zero CoHex density^+ in the minor groove, so the corresponding panoramic density is not shown. . . ... [TO] 

6 Simulated and experimental ASAXS profiles for 25 bp DNA in (a) 100 mM RbCl and (b) 10 mM SrCl 2 solu¬ 

tions. Experimental data (ITOllTOl) are shown as black dots. This hgure shows simulation results using the full 
3D cDFT model (blue lines), the 3D cDFT model without ion-water interactions (red lines), the full ID cDFT 
model (blue dots), and the NLPB model (green lines and dots). . dn 

7 Simulated and experimental ASAXS profiles for 25 bp DNA in 0.5 mM CoHexCla DNA solutions. Experi¬ 
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line. ID cDET results are shown with a thin black line and ID NLPB results with a red line. The 3D cDET 
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